\documentclass{ctexart}
\usepackage{graphicx}
\begin{document}

\section{\label{sec1}计算地月运输轨道}

在短时间内，近似的认为月球轨道为圆形。
计算地月轨道所使用的参考系称作$S_1$：以地月系统的质心G为坐标原点，与月球的公转一同转动的运动坐标系$x_1 y_1 z_1$，其中$x_1 y_1$平面是月球公转平面，$x_1$轴为地月连线，如图：

设$m_1$、$m_2$分别代表地球和月球，m代表小物体，Ω为月球公转角速度，$r_1$、$r_2$、r的意义如图，$\pi_1$、$\pi_2$分别为$\pi_1=\frac{m_1}{m_1+m_2}$, $\pi_2=\frac{m_2}{m_1+m_2}$。则卫星运动的限制性三体问题的标量方程为：
\begin{equation}
\left\{ \begin{array}{l}
\ddot{x}=2\Omega\dot{y}+\Omega^2x-\frac{Gm_1}{r_1^3}(x+\pi_2 r_{12})-\frac{Gm_2}{r_2^3}(x-\pi_1 r_{12})\\
\ddot{y}=-2\Omega \dot{x}+\Omega^2 y-\frac{Gm_1}{r_1^3}y-\frac{Gm_2}{r_2^3}y\\
\ddot{z}=-\frac{Gm_1}{r_1^3}z-\frac{Gm_2}{r_2^3}z
\end{array} \right.
\end{equation}

验证能量守恒用的方程为
\begin{equation}
\frac{E}{m}=\frac{1}{2}v^2+\frac{1}{2}\Omega^2r^2-\frac{Gm_1}{r_1}-\frac{Gm_2}{r_2}
\end{equation}


\section{\label{sec2}参考系之间的转换}

\includegraphics{graph1.png}

\includegraphics{graph2.png}

\includegraphics{graph3.png}

\includegraphics{graph4.png}

\includegraphics{graph5.png}

\includegraphics{graph7.png}

\includegraphics{graph8.png}

$S_1$：以地月系统的质心G为坐标原点，与月球的公转一同转动。其中$x_1 y_1$平面是月球公转平面，$x_1$轴为地月连线。

$S_2$：原点为地球的地心，$x_2 y_2 z_2$方向始终同$x_1 y_1 z_1$一致。
因此有
\begin{equation}
\left\{ \begin{array}{l}
x_1=x_2-r_{12}*\pi_2
\\y_1=y_2
\\z_1=z_2
\\v_{x1}=v_{x2}
\\v_{y1}=v_{y2}
\\v_{z1}=v_{z2}
\end{array} \right.
\end{equation}

$S_3$：原点为地球的地心，$z_3$与$z_2$始终一致，$x_3$与$y_3$为平动参考系，并不随着月球的公转而旋转。因此有
\begin{equation}
\left\{ \begin{array}{l}
x_2=x_3\cos\gamma-y_3\sin\gamma
\\y_2=x_3\sin\gamma+y_3\cos\gamma
\\z_2=z_3
\\v_{x2}=v_{x3}\cos\gamma-v_{y3}\sin\gamma+\Omega*R_{earth}*\sin\gamma
\\v_{y2}=v_{x3}\sin\gamma+v_{y3}\cos\gamma-\Omega*R_{earth}*\cos\gamma
\\v_{z2}=v_{z3}
\end{array} \right.
\end{equation}
(其中$\gamma=\gamma_0-\Omega*t$，$\Omega$为月球公转角速度,$R_{earth}$为地球半径。)

$S_4$：原点在地心，$y_4$与$y_3$相同，$z_4$轴为地球的地轴。$\alpha$为$z_3$与$z_4$的夹角，即黄赤交角。则
\begin{equation}
\left\{ \begin{array}{l}
x_3=x_4\cos\alpha-z_4\sin\alpha
\\y_3=y_4
\\z_3=x_4\sin\alpha+z_4\cos\alpha
\\v_{x3}=v_{x4}\cos\alpha-v_{z4}\sin\alpha
\\v_{y3}=v_{y4}
\\v_{z3}=v_{x4}\sin\alpha+v_{z4}\cos\alpha
\end{array} \right.
\end{equation}

$S_5$：原点在地心，$z_5$与$z_4$相同，$x_5 y_5$则跟随地球的自转一同旋转。设地球上P点的位置参数$\phi\theta$如图。则有
\begin{equation}
\left\{ \begin{array}{l}
x_4=x_5
\\y_4=y_5
\\z_4=z_5
\\v_{x4}=v_{x5}-\omega_{earth}R_{earth}\cos\theta\sin\phi
\\v_{y4}=v_{y5}+\omega_{earth}R_{earth}\cos\theta\cos\phi
\\v_{z4}=v_{z5}
\end{array} \right.
\end{equation}

设发射速度为u，发射的方位角$\xi\eta$如图，则有
\begin{equation}
\left\{ \begin{array}{l}
x_5=R_{earth}\cos\theta\cos\phi
\\y_5=R_{earth}\cos\theta\sin\phi
\\z_5=R_{earth}\sin\theta
\\v_{x5}=u\sin\eta\cos\theta\cos\phi-u\cos\eta\cos\xi\sin\phi-u\cos\eta\sin\xi\cos\phi
\\v_{y5}=u\sin\eta\cos\theta\sin\phi+u\cos\eta\cos\xi\cos\phi-u\cos\eta\sin\xi\sin\phi
\\v_{z5}=u\sin\eta\sin\theta+u\cos\eta\sin\xi\cos\theta
\end{array} \right.
\end{equation}

最终完成从$(u,\xi,\eta,\phi,\theta,\gamma)$到$(x_1,v_{x1},y_1,v_{y1},z_1,v_{z1})$的转化。

\section{\label{sec3}空气阻力}

空气阻力公式为$F=\frac 12 C_d \rho A v^2$。在高超音速的情况下，Cd需要做详细计算。最终的计算结果发现，只要形状设计的足够好，例如一个$\theta=10^{\circ}$的尖锥形状飞行器，Cd可以做到0.1那么小。当面积$A=1m^2$，质量$M=10t$时，程序算得要飞出大气层，初速度需要9.8km/s，而飞出大气层时，速度只有初速度的$35\%$。

（以下内容参考《高超声速飞行器空气动力学》黄志澄 编著）：总阻力系数分为以下四个部分：无粘阻力系数、摩擦阻力系数、粘性干扰阻力系数和底部阻力系数（见P113）。

\subsection{\label{无粘阻力系数}无粘阻力系数}

在P112上写道，根据实验数据拟合出来的机身K值为
\begin{equation}
K=2.38+0.03792\theta-0.002521\theta^2+0.00004583\theta^3+2.917\theta^4
\end{equation}
其中K值的含义为$C_p=K\sin^2\theta$，而$C_P$的含义则是$\frac {F_p} A=\frac12 \rho C_p v^2$. $F_p$就是无粘阻力，A是截面积。其中$\theta$是形状有关的角度，可见飞行器越尖，这部分阻力越小。为了简单，可以假设我们的飞行器是尖锥形状的。

\subsection{\label{摩擦阻力系数}摩擦阻力系数}

根据可压缩边界层的积分方法，在湍流边界层时，尖锥的$C_f$与平板的$C_f$相差不大(见P207)：
\begin{equation}
\frac{C_f}{C_{f-board}}=1.086~1.176
\end{equation}
其中$C_f$为摩擦阻力系数，定义为$\frac {F_f} S=\frac12 \rho C_p v^2$，其中$F_f$为摩擦阻力，S是侧面积。

而平板的$C_f$则可以根据《空气动力学》（陆志良 等编著）P125页的经验公式计算：

当$Re>10^7$时，
\begin{equation}
C_f=\frac {0.455} {{(lg Re)}^{2.58}}
\end{equation}
其中Re为雷诺数：$Re=\frac {\rho v l}{\mu}$，l为平板长度。
要注意，虽然$C_f$与Re为负相关，但是在算受力的时候实际上就变成正相关了。最后计算时要注意不要在$rho->0$的时候弄出来个零除以零就好。。

\subsection{\label{粘性干扰阻力系数}粘性干扰阻力系数}

没看懂这部分怎么搞。。是不是已经包含到上面那部分里面去了？、、

\subsection{\label{底部阻力系数}底部阻力系数}

在P113页上写道，底部阻力系数随Ma数的上升迅速下降，在我们的情况中，这一项相对于前几种阻力来说可以忽略不计。

\section{\label{sec4}计算最优轨道}


采用的算法为遍历搜索法，先确定大范围，再逐步缩小搜索范围。一开始时为了避免错过解，先把月球半径扩大了来搜索解。

在计算时，我们认为$\phi$与$\gamma$并不那么独立，因此事先规定$\phi=-\frac{\pi}{2}$。$\theta$代表纬度，因为发射高度对速度影响较大，因此事先规定好发射地点在青藏高原。因此事先规定$\theta=0.138879\pi$，即青藏高原的纬度。

最终计算得的结果是，$u_0$最小为9800m/s，$\eta_0=0.289063\pi$，$\xi=-0.214844\pi$，$\phi=-0.500000\pi$，$\theta=0.138879\pi$，$\gamma=0.347656\pi$，单位全部为弧度。

\section{\label{sec5}电磁加速}


利用公式
$F=\nabla(m\cdot B)=m \frac{\partial}{\partial x}B$
其中m为抛射体的磁矩,$m=NIA$。认为抛射体磁铁可以等效为有限长通电螺线管，则它自己产生的磁场可以近似为为$B_m=\mu_0\frac N l I \frac l {\sqrt{l^2+{(2r)}^2}}$，目前永磁体能产生的最大磁场约为$B_m=1.4T$。因此$m=B_m \sqrt{l^2+{(2r)}^2} A/\mu_0$。其他参数为l=10m，$A=1m^2$，r为半径。
我们认为m为一常数，因此加速度公式就可以简化为$a=\frac{F}{M}=\frac{m}{M}\frac{\partial}{\partial x}B$

根据电磁学知识，我们知道环形线圈在轴线上产生的磁场为$B_x=\frac{\mu_0 R_{coil}^2 I}{2(R_{coil}^2+x^2)^{3/2}}$，其中：Rcoil为加速线圈的半径，可以设为Rcoil=1m；x为x方向上离线圈的距离。

然后I(t)和B(x,t)使用LRC放电电路来计算，见第6部分。

设线圈间距为lcoil。设抛射物中心的x坐标为h。则有$h(T)=\int\limits_0^T v(t)dt+h_0$，而$v(t)=v_0+\int\limits_0^t a(\tau)d\tau$。
因抛射体每一个线圈的位置不同，因此需要把整个螺线管拆成一个一个的线圈来处理。设dm=m/N，dx*N=l，所以dm=m/N=m/l*dx。
因此$a(t)=\int\limits_{h(t)-l/2}^{h(t)+l/2}\frac{dm}{M} \frac{\partial}{\partial x}B(x)=\frac{m}{Ml}\int\limits_{h(t)-l/2}^{h(t)+l/2}\frac{\partial}{\partial x}B(x)dx =\frac{m}{Ml}(B(h(t)+l/2)-B(h(t)-l/2))$。
用数值方法可以解出上面的a(t)及h(t)。

\section{\label{sec6}电容放电}

认为电容放电是个LRC电路。有关LRC电路放电的电流，由下图的公式所确定：

\includegraphics{currency.png}

要得到一个确定的电流曲线，我们需要知道L、C、R、U0四个指标。我们现在想要做的，是找出这四个指标可以达到的一个组合，使得极大放电电流最大，从而得到最大的磁场。为此，我们需要反推：现有技术最大的放电电流能产生50T的磁场。而通过考察极大电流与这四个变量的关系可知，L、R越小，C、U0越大，则极大电流值越大。通过这个逻辑，我们可以定出一组数据以产生相应的磁场：L=0.1mH, C=0.1F, R=5m$\Omega$, U0=5MV，此时最大磁场约为42T。

\includegraphics{curve_of_B.png}

\includegraphics{curve_of_a.png}

关于如何同步：我们的想法是让抛射体的前端恰好达到一个线圈时，这个线圈应当处于其电流最大值处，因此$h_0=-l/2-T/4*v0$（T为振荡电路的周期）。另外还得注意，电流波形里sin中的那个w必须很大，否则当抛射体穿过线圈之后那个磁场就成了起阻碍作用的磁场了。

\section{\label{sec7}结论}
根据$v^2=2as$，我们要让出射速度为10km/s，还要让轨道控制在5000m左右，那么需要的平均加速度是1000g。上面的数据貌似还没达到那么大。最关键的是，内禀矫顽力的限制使得外磁场也不可能那么大还让磁体保持原状。

\section{\label{sec8}抛射物为超导体}
若抛射物为超导体，则需要根据迈斯纳效应算受力，而不再是之前的公式。迈斯纳效应受力的计算为，先算出为了保持内部磁场为零所需要的电流，然后这个电流在磁场中的受力即迈斯纳效应导致的受力。这个力与xy方向的B有关，而不是Bz了。因为磁铁也可以等效为一个环形电流，因此实际上超导体受力和磁铁是一样的，只是需要注意第二类超导体并不具有完全的抗磁性，前面需要乘以一个小于1的系数。

\section{能量利用效率}
若按照5000m轨道上每3m放置一个线圈的话，储存的总能量应为
\begin{equation}
E_{coiltot}=\frac12 C U_0^2*N=8.3*10^{12}J
\end{equation}
而抛射体出射时的动能为
\begin{equation}
E_k=\frac12 M v^2=5*10^{11}J
\end{equation}
能量利用率6\%。。。

\section{抛射体为互感线圈时的发热问题}
按照内部感应出的电流为兆安量级来算。铜在20℃时的电阻率为$8*10^{-8}\Omega m$，铜线圈直径按照5cm算，10m长得话，电阻大约为$0.1m\Omega$。这样发的总热量为$I^2Rt=10^{8}J$。按照密度$8.9*10^3kg/m^3$来算，质量大约在700kg。铜比热容为$0.39×10^3J/(kg·℃)$，算得升温400℃左右。
\begin{equation}
\Delta T=\frac{I^2t\rho_{elec}}{c\rho_{mass}S^2}
\end{equation}
铜的熔点在1000℃左右。这个导线有点太粗了。。但还是有可能的。。。

\section{\label{sec9}各个文件的含义作用}
magnet4.m：给定v0、L、C、R、U0、l、lcoil之后，算出I(t)以及B(x,t)，然后根据磁场位形计算出在一个线圈的作用下抛射体的加速度随t或者h如何变化。

drag.m：给定初始速度以及仰角，算出由于空气阻力的作用，在抛射体射出大气层后的速度、仰角以及速度变为原来的多少倍。

calculate.m：给定物体的初始位置和初始速度，计算出它在地球、月球引力下的轨迹。

coordinate.m：参考系转换用的。给定初始速率和一堆角度，算出calculate.m所需要的初始位置及初始速度。

main2.m：遍历各个角度，用来搜寻最小发射速率。

currency.nb：算出LRC电路的电流曲线。

LCRU0.m、accelerate3.m、draw.m：貌似已经发挥完了作用了，现在已经没用了。









\end{document} 